Mathematical Modeling and Numerical Simulation for the Outbreak of COVID-19 Involving Loss of Immunity and Quarantined Class

For the analysis of the recent deadly pandemic Sars-Cov-2, we constructed the mathematical model containing the whole population, partitioned into five different compartments, represented by SEIQR model. This current model especially contains the quarantined class and the factor of loss of immunity. Further, we discussed the stability of the SEIQR model (constructed on the basis of system of coupled differential equations). The basic reproduction that indicates the behavior of the disease is also estimated by the use of next-generation matrix method. Numerical simulation of this model is provided, the results are analyzed by theoretically strong numerical methods, and computationally known tool MATLAB Simulink is also used for visualization of the results. Validation of results by Simulink software and numerical methods shows that our model and adopted methodology are appropriate and accurate and could be used for further predictions on COVID-19. Our results suggest that the isolation of the active cases and strong immunization of patients or individuals play a major role to fight against the deadly Sars-Cov-2.


Introduction
At present, emergence of deadly respiratory virus COVID-19 throughout 209 countries of the world is a major global concern. Primarily, this virus was known as Severe Acute Respiratory Syndrome Coronavirus-2 (Sars-Cov-2), which initially came out in Wuhan city of China. COVID-19, which originated at the end of December 19 in Wuhan, China, is considered to be the third epidemic of CoV, and it is holding almost the same symptoms like Sars-CoV. And this disease is found to be more deadly than the coronavirus Sars.
Since its outbreak, this virus has caused enormous deaths up to 4,182,831 and has infected 195,345,791 people worldwide. The main cause of the outspread of this virus is the contact of infected person with healthy individuals because it was found in study that this infection is usually caused by transmission of globules through coughing or sneezing. These globules can stay in the air for a long time and can cause infection to others. However, it is quite challenging for the scientists to investigate the preventive measures under which the spread of this virus can be controlled and to produce a vaccination to fight against this virus. Until now, almost 13 vaccinations against this disease exist and scientists are still devoting their attentions to produce a strong vaccination against this virus. A huge amount of research has been carried out to look over the conditions and circumstances by which this deadly virus can be controlled. Scientists find out COVID-19 to be one of the crucial outbreaks that attack the respiratory system [1]. One of the main reasons behind the outspread of COVID-19 is due to the transmission of germs through respiratory globules among humans, and this virus is considered to be the vector transmission. the WHO (World Health Organization) warned that if control measures are not implemented in time, then the outbreak of the coronavirus can spread more rapidly [2].
Qualitative analysis using the concept of basic reproduction ratio estimated by next-generation matrix and stability theory of differential equations to examine the outburst of COVID-19 was studied by Olumuyiwa et al. [3]; basically, they estimated the epidemic behavior in Pakistan by using some previous statistical data from Pakistan and they applied numerical approaches Nsfd and Ode45 to figure out the study.
Researchers are adopting different tactics to investigate the growing behavior of the coronavirus and carrying out their studies to find out the ways to bring COVID-19 to an end. One of the methods of analyzing the behavior of diseases is compartmental modelling that is applicable to the mathematical models involving influential diseases. In this approach, the population is allocated to various compartments with particular labels, for example, susceptible (S) and infected (I). Mathematical models are quite helpful in investigating the behavior of viruses, diseases, or infections and help out to conclude under what circumstances the outbreaking virus can be eradicated or continued in citizens and are often useful in estimating the duration of an outbreak. These models are usually considered to take form of differential equations.
Annas et al. constructed the mathematical model to study the influential behavior of coronavirus [4,5]; they also discussed the stability of their proposed model. Suleman et al., Ul Rahman et al., and Iqbal and Karaca numerically investigated the fractional model of HIV [6,7,8]. The National Health Organization instructed to keep away from infectious person or animals either with fever or any respiratory problems and also directed to wear surgical masks and directed continual use of hand sanitizers in public places for self-safety from infection. Social distancing or less crowded places can reduce the risk of spreading of coronavirus because it is more likely to spread in compact places. Batista constructed the SIR model to guess the finishing measurements of coronavirus [9]. Moreover, the basic reproductive number ðR 0 Þ is very effective in approximating the transmission rate of an infection; i.e., it is also useful in estimating the ratio of occupants required to be immunized in order to wipe out the infection.
Basic reproductive number is estimated basically from the mathematical model, and in most epidemic models, the spreading rate of infection or disease-free equilibrium point is said to be endemic or stable if R 0 lies between 0 and 1, otherwise it will be epidemic or unstable. Substantially, the larger the values of R 0 , the harder it is to take control of the outbreak of the disease. Zhao et al. proposed the test approximation of the basic reproductive number ðR 0 Þ for the breakout of COVID-19 during its early stages [10,11]. However, it is not that easy to find out the explicit expression for the basic reproductive ratio so more advanced approaches, e.g., next-generation matrix, are required to estimate the reproductive number [12]. The study of basic reproduction ratio ðR 0 Þ for dissemination of epidemical infection was presented by Van den Driessche and Watmough, and they also bring forward the stability and unstability of infection free equilibrium based on reproductive ratio [13,14].
The use of epidemic computational simulation models play a key role in estimating the transferal parameters and in guessing the influential behavior of the infection. These models are proven to be useful in depicting the growth rate or decay rate of viruses with the passage of time and are quite helpful in providing the control measures that can be adapted to reduce the spread of the disease. A lot of results and applications carried out from the simulation models of COVID-19 are being accepted and published. Abdulrahman presented the computer software Simulink program to track the contagious virus COVID-19. He used SEIRD and SIR epidemic models in the form of algebraic-differential equations and simulated the results by using Simulink approach [15].
One way of examining the influential growth of this deadly outbursts is to use the computer simulation block models. In various publications, different numerical as well as analytical approaches are taken into account to carry out the outputs of SIR, SEIR, and SEIRD models. A SIR model involving the immigration rate was proposed by Ud Din et al. They simulated their outputs by applying the approach of Nsfd and investigated the changing behavior of epidemic due to immigration factor between different classes [16].
Numerical approaches are also useful in analyzing the changing behavior of the outburst diseases and can provide some better precautionary measures to prevent the diseases transmission. In literature, average researchers have used meshless as well as numerical methods like FDM, FVM, FEM, Euler, and Runge-Kutta. But the classical numerical methods like RK and Euler are less computational, efficient, and easily applicable, as meshless methods are useful in simulating physical phenomena including biological as well as engineering-related problems, as if analyzed the SEIR model of disease by applying meshless methods like EFBMM and OSBMM [17].
Ahmed et al. adapted some numerical techniques, including Rk2, Rk4, and Euler's method as well as simulation process to study some mathematical models developed for COVID-19. Further, they estimated the epidemic size for Iraq and Turkey with the help of a logistic model [18,19]. Maier and Brockmann presented the growth rate of COVID cases in contrast with the initial rate of suspected cases [20,21]. Baba et al. numerically investigated the effects of the lockdown by constructing the model based on five equations and also discussed the equilibria and its stability [22,23]. Zha et al. proposed the fuzzy-based approach in order to lessen the outbreak of the novel Sars-Cov-2 [24].
Bassetti et al. studied the contaminating behavior of COVID-19, and further, they investigated what obstacles we can face due to this virus [25,26]. Cao [34]. Khoshnaw et al. pointed towards the significance of using computational simulation tools to forecast the behavior of infections. They bring forward the idea of sensitivity analysis to test the sensitivity of the models formulated on the basis of differential equations [35,36]. Ul Rahman et al. proposed the model for paint industry effluent and simulated the results by using MATLAB Simulink tool [37].
The short-term analysis of Sars in three huge cities of India was investigated by Mandal et al.; also, they advanced the method to control out the basic reproduction ratio [38,39]. Liu  Egbetade et al. analyzed the SIR model of the infections and also discussed the existence of equilibria and the reaction of disease on the basis of R 0 [45,46]. Ul Rahman et al. used the concept of numerical simulation to model some industrial problems and analyzed the models by using numerical techniques and Simulink [47,48]. Abdulrahman, Iqbal and Wu, and Rahim et al. used simulation programs to analyze the mathematical models constructed to study COVID-19 and biological models [49,50,51,52].
From the above studies, it has been concluded that the reason behind the global outspread of COVID-19 is dense places and migration of infected people and loss of immunity in people to fight against diseases. Therefore, to control the outspread of CoV, it is necessary to quarantine the infected person and moreover, strong immunization of patients must be necessitated to get rid from this infection.
(i) In this study, we will estimate the control of some parameters including immunity parameter and isolation class that will be beneficial to provide some controlling features for the above-mentioned disease. This research work is constructed by the SEIQR model based on nonlinear differential equations, and the model is analyzed by the simulation block models and by using some numerical approaches. Furthermore, the stability of the reproduction ratio (R 0 ) and the existence of disease-free equilibria are also studied. The idea of next-generation matrix is utilized to find out the reproduction ratio (ii) The summary of our work is as follows: we applied Euler's method, Rk4, and Ode45 to obtain the results for our differential-equation epidemic model (iii) Another aim is to estimate the local stability of reproductive number For this purpose, we divided the paper in three sectors. The first section is for the mathematical model; in the second section, the stability analysis and the existence of equilibrium point are discussed; and in the last section, the numerical and graphical outcomes are represented.
The flow diagram of the SEIQR model is presented below. In the flow diagram (Figure 1), the transmission of population with different transition rates, among the five compartments, is shown. In Figure 1, the parameter Z is representing the rate of population joining the susceptible class and further, the individuals from the susceptible community are moving to the infected as well as exposed community with β force of infection. The parameter πE indicates the switching rate of individuals from the exposed to the infected class. The isolated class contains the individuals from the exposed and infected class, with γ and σ joining rates, respectively. Further, recovered individuals with lack of immunity leave the recovered class with α rate and move to the susceptible class.

Formulation of Mathematical Model
The SIR model is the simplest fundamental model constituting three compartments of complete population. And it was first used in 1916 and then in 1927 to estimate the behavior of diseases and viruses. Other models like SEIRD, SIRD, and SEIR are the extensions of this basic model.
The present model is the SEIQR model. As the reason behind the upsurge COVID-19 is the lack of immunity and the contact of infected patients with other healthy individuals, therefore, in the said work, the SEIQR stochastic mathematical model is formulated involving the immunity parameter and the quarantined community.
For this, the partition of whole population is placed in five different compartments named as (i) susceptible (S), (ii) exposed (E), (iii) infected (I), (iv) quarantined (Q), and (v) recovered (R).
The susceptible class contains those individuals having mild symptoms and who are at risk of getting infectious Table 1. The individuals who are needed to be quarantined are indulged in isolated class, and people who caught the disease are present in the infected class, whereas the recovered compartment contains those individuals who are either dead or have recovered from infection or the people who are still at risk of getting infectious again due to less immunity.

SEIQR Mathematical
Model. The SEIQR model is given by the following nonlinear 1st-order differential equations, and the description of each compartment is given Table 1.

Computational and Mathematical Methods in Medicine
Each equation is describing the transmission behavior of individuals in the respective compartments. By this transmission, number of individuals can vary in each of the five compartments. The description of the transition rates in each cell is given in Table 2.
Let us define the initial conditions to be Sð0Þ The precise definition of the compartments used in the formulation of the model can be seen from Table 1

For Positivity and Boundedness of Solution.
For the positiveness of the solution and bounded solution of the above system, it is necessary that the solution maintains nonnegativity for all t ≥ 0.
And it was found that The interpretation of the parameters used in the formulation of the model can be seen from Table 2 given below.

Compartments
Brief definition S Susceptible community (susceptible to disease) E Exposed community (those people who come in contact with the virus) I Infected community (when someone is exposed to the disease and having 70% symptoms) Q Quarantined community R Recovered community To determine R 0 , we used the next-generation matrix method and obtained since R 0 is defined by the relation R 0 equals the eigenvalue of FV −1 , where F includes the terms with secondary infectious disease and V includes the terms other than the secondary infectious disease.
And the stability of system can be decided from R 0 , i.e.,   The stability discussed here is basically local stability.

Simulink Block Model.
In order to speculate and trace the eruption of Sars-Cov-2, we proffered the computer-based simulation scheme in this considered article. For this objective, the respective differential equations are simulated with the help of blocks available in the library browser of Simulink-MATLAB. Simulink tool is easy to use for the purpose of predicting behavior of any natural phenomenon or system. The above considered mathematical model is rooted in Simulink with the help of block diagram given in the appendix.

Numerical Simulation and Results
For the output of the present model, we implemented two numerical methods Euler's method and Rk4. Also, we   The dynamical behavior of population in all five classes is shown graphically by taking different time periods. Further, the comparison of three methods Ode45/Simulink, Rk4, and Euler method is shown in graphs.
The changing behavior of population in various compartments for the stable infection, i.e., for endemic case R 0 < 1.
Dynamics of population over 365 days, in various compartments for unstable infection, i.e., for pandemic case R 0 < 1.
As no proper vaccination for the disease is discovered yet, so the entire individuals are at risk of getting contaminate by this infection. That is why the whole population is often put into the susceptible class. From the susceptible class, the entities of that class can get infected by making contact with an infectious person and then, they can join the infected class as well as the exposed class. The differential behavior of the susceptible entities is shown in Figure 2(a) for the case of endemic disease. In Figure 2(a), the behavior is depicted over the period of 20 days as one can see that when we have an endemic case,    7 Computational and Mathematical Methods in Medicine so the susceptible class has a constant rate of variation as people are at less risk of getting engaged to the disease and there will be less movement of individuals from this class to other compartments. Moreover, in Figure 2(b), the graph is depicting the changing behavior of the population in the susceptible class over the period of 40 days and in Figures 3(a) and 3(b), the variation in susceptible individuals is shown over the period of 60 and 80 days, respectively. Another aspect of the constant variation of individuals in the susceptible class is that when we have the population with strong immunity system, i.e., less people will get infected to the disease. The abovementioned graphs are drawn for the strong immunity parameter. And all these results are carried out with the help of Ode45, Rk4, and Euler's method.
In Figures 4(a) and 4(b), the transition behavior of the population having less immunity is discussed, since for less     Computational and Mathematical Methods in Medicine immunization factor, there is higher probability of getting involved with the infection. As one can see from both graphs, individuals from the susceptible cell are becoming part of infected as well as exposed cells. The decreasing behavior is depicting the tendency of people getting infectious and exposed. The variation of entities in the exposed cell considering the case of endemic virus is depicted in Figures 5(a), 5(b), 6(a), and 6(b). One can see that the graph of this class is illustrating the sharp reduction throughout different periods of time. The reason behind this decrement is the immigration of individuals from the exposed class to the infected class. As the exposed individuals have almost 60% signs of the disease, so they are almost considered to be infected and therefore, the graphs are showing the sharp rate of variation from the exposed cell to the infected one.    The evolution of individuals in the infected community is portrayed in Figures 7 and 8. The transition behavior of individuals in the infected class is depicted in Figure 7(a) over the period of 20 days; the initial increment is due to the joining of the infected cell of exposed people, and in Figure 7(b), the graph is drawn for the time period of 40 days; the early increasing behavior is due to the increase of population in the infected cell, and the decreasing behavior later is due to the movement of patients from the infected to the quarantined cell. In Figures 8(a) and 8(b), the dynamical behavior of the population is illustrated over the period of 60 and 80 days, respectively. The diminution after some time is due to the immigration of patients to the isolated cell.
As we are considering the isolation of patients in order to prevent the further spread of the disease, so the infected patients are needed to be quarantined for the reduction of   the spread. The dynamical behavior of the quarantined population is shown in Figures 9 and 10 for different time periods. From these graphs, it can be seen that people are getting involved in the recovered cell with sharp rate if they are quarantined in time. The variation of recovered population is drafted in Figures 11 and 12; it can be seen that the graphs are depicting an increasing behavior; the reason behind this sharp recovery is due to the isolation of infected individuals and the strong immunization of the patients. The initial decrement is due to the less immunity factors and the deaths of patients either by disease or other reasons.
The dynamics of population in the susceptible class is shown in Figures 13(a) and 13(b); in Figure 13(a); the behav-ior is shown for the strong immunity parameter; on the other hand, the graphical results are depicted by considering weak immunity among individuals. In Figure 14(a), transition behavior in the exposed class is depicted over 365 days; the decreasing effect is due to the transition rate of individuals from the exposed class to the infected class, in which in Figure 14(b), the dynamics of population in the isolated class displays the initial increment and then decrement of population; this fluctuation is due to the joining and motion of individuals in this class from other classes. In Figure 14(c), the dynamics of individuals shows the rise in population of the isolated class and the later declining behavior shows the flow population from the isolated class   Figures 15(a) and 15(b), the weak and strong immunity among individuals is considered, respectively; Figure 15(a) shows a slow rate of recovery due to loss of immunity; and Figure 15(b) shows a sharp rate of recovery due to strong immunity. Now, the illustration for the transmission dynamics of population in each compartment is considered for the case of epidemic disease. In Figures 16(a) and 16(b), the rate of change of population in the susceptible class is illustrated for 20 and 40 days, respectively, by taking into account the epidemic behavior of disease. We can see the suppressed behavior of population in the S class that is due to the fact that people from this class can catch the infection and can move to other compartments. In the endemic case, we have a higher rate of transition from one cell to another, as compared to the epidemic case. The same reaction of population   is shown in Figures 17(a) and 17(b) considering 60 days and 80 days, respectively.
In Figures 18 and 19, transmission dynamics of population in the exposed cell is given. The prime rise is by the increase in population of this cell, as from the susceptible class, more individuals are engaging in this class and the later decrease is due to the fact of epidemic behavior, as disease is not in control so exposed individuals are getting involved in infected cell. And the reason behind the inflation of population in the infected class is the excessive amount of indulgence of individuals between the exposed and infected classes due to the epidemic behavior as one can see in Figures 20 and 21. Figures 22 and 23 show the dynamics of population in the isolated compartment over 20, 40, 60, and 80 days, respectively. The recovered rate of the population is presented in Figures 24 and 25; it can be seen from the graphs that in the case of epidemic virus, we have slow-going recovery of infected individuals as compared to those in the case of endemic behavior. Similarly, the dynamical behavior of the virus is endemic and stable so less people will catch the infection (ii) for  Figure 19: Dynamical behavior of population in the exposed compartment.

18
Computational and Mathematical Methods in Medicine the virus is epidemic and unstable so, more people will catch the infection The changing behavior of population in various compartments for unstable infection, i.e., for pandemic case R 0 > 1.
Dynamics of population over 365 days, in various compartments for unstable infection, i.e., for pandemic case R 0 > 1.

Conclusion
In the said study, we concluded that isolation of the infected person and strong-immunization of the individuals can be taken as precautionary measures to control the virus. By following the Sops implemented by the government, we can prevent the further outbreak of this disease. For the unstable case, it can be seen that there is sharp variation of entities among all five classes. And for the lower immunity factor, individuals are still at risk of getting infectious and there are fewer chances of recovery. And for the stable case, it is apparent from graphical results that less people are exposed to infection. As in the susceptible class, we have constant behavior of population with zero loss-immunity parameter. But if we include the lack of immunity parameter, we can see that more people are joining the exposed and isolation classes but the situation is in control since we have endemic behavior here. In Section 3, it is discussed in detail. The comparison of the three methods can also be seen from the graphs. It can be seen that Ode45 is more accurate than Euler's method (1st-order RK method) and Rk4 as it is built in MATLAB and it implements Rk4 or Rk5 to solve the equations.
Validation of results by Simulink software and numerical methods shows that our model and adopted methodology are appropriate and accurate and could be used for further predictions on COVID-19.
Simulink software is helpful to analyze the performance of any system or reaction of any natural phenomenon. Simulink tool is quite handy for the prediction and forecasting of any disease or infection. This tool is based on block diagrams and is beneficial in modeling and investigating the dynamical behavior of various phenomena. Engineers and scientists use Simulink software to test their proposed designs and models before finalizing their results. A mathematical model (constructed with the help of block diagrams) is required to start up the simulation program.
It is available in the main page of MATLAB. One can click on the Simulink icon, and it opens up with the Simulink start page. There, it appears with Simulink library browser constituting a huge library of blocks like math operations, sinks, sources, and commonly used blocks. The above mathematical model presented in Figure 29 is simulated by the formulation of block models based on integrators, gains, add blocks, product, and scope.
The product, integrator, sum, gain, and scope are available in the commonly used block library but they have different libraries too. The add and gain block is also available in the library of math operations; display and scope can be found in sinks. The add block is used to sum the inputs, the product block is used the multiply the inputs, and scope represents the results graphically. And integrator is accessible from continuous block library. The integrator block integrates the input signal with respect to time.
The outputs for the presented formulated model was simulated with the help of computer simulation software which solves the models with the help of block diagrams. The Simulink model is given by the following figure.

Data Availability
The numerical data used to support the findings of this study are included within the article.

Additional Points
Highlights. The SEIRQ model is modified to study the epidemic Sars-Cov-2. The problem is analyzed by MATLAB Simulink software and numerical methods (Euler method) and ðRk4Þ. Evaluation of basic reproductive number is done by the next-generation matrix method. Local stability of infection-free equilibria is discussed. Numerical simulation is performed for each class considered in the modified SEIQR model. Results obtained for each class considering epidemic and endemic behavior of disease are represented graphically